%DG_1D Burger u(x,0) = 0.5+sin x, simplify
clc ; clear ; close all;
h = [2*pi/20,2*pi/40,2*pi/80,2*pi/160];
errors = zeros(4,1);
k=2;
for d = 1 : 4
    X = (0:h(d):2*pi);
    x_half = 0.5*(X(1:end-1)+X(2:end));
    cfl = 0.2;
    N = length(X)-1;
    u0 = zeros(N+1,1);
    u_true = zeros(N,1);
    %生成基函数
    [b_fun,db_fun] = pro_basis(k,X);
    for i = 1 : N+1
        u0(i) = 0.5+sin(X(i));
    end
    %求初始系数
    A0_coef = zeros(k+1,N);
    A = zeros(k+1,k+1);
    % for i = 1 : N
    %     x_x = linspace(x(i),x(i+1),k+2);
    %     x_x = x_x(2:end);
    %     b = zeros(k+1,1);
    %     for l = 1:k+1
    %         Pl = b_fun{i,l};
    %         for j = 1 : k+1
    %             A(j,l) = Pl(x_x(j));
    %         end
    %         b(l) = 0.5+sin(x_x(l));
    %     end
    %     A0_coef(:,i) = A\b;
    % end
    for i = 1 : N
        M = zeros(k+1,k+1);
        b = zeros(k+1,1);
        for m = 1 : k+1
            Pm = b_fun{i,m};
            for l = 1 :k+1
                Pl = b_fun{i,l};
                PmPl = @(x) Pm(x).*Pl(x);
                M(m,l) = quadgk(PmPl,X(i),X(i+1));
            end
            uPm = @(x) (sin(x)).*Pm(x);
            b(m) = quadgk(uPm,X(i),X(i+1));
        end
        A0_coef(:,i) = b;
    end
    T = 0;
    u0 = zeros(N,1);
    for i = 1 : N
        u0(i) = sin(x_half(i));
    end
    alpha = max(abs(u0));
    dt = (0.19*h(d))/alpha;
    while T <0.2
        T = T+dt;
        [New_u0,New_A0_coef] = Euler_forward(u0,A0_coef,X,x_half,k,N,h(d),dt,b_fun,db_fun);
        [u1,A1_coef] = Euler_forward(New_u0,New_A0_coef,X,x_half,k,N,h(d),dt,b_fun,db_fun);
        u2 = 0.75*u0+0.25*u1;
        A2_coef = 0.75*A0_coef+0.25*A1_coef;
        [u3,A3_coef] = Euler_forward(u2,A2_coef,X,x_half,k,N,h(d),dt,b_fun,db_fun);
        u0 = (1/3)*u0+(2/3)*u3;
        A0_coef = (1/3)*A0_coef+(2/3)*A3_coef;
        alpha = max(abs(u0));
        dt = (0.19*h(d))/alpha;
        % u0 = 0.5*u0+0.5*u1;
        % A0_coef = 0.5*A0_coef+0.5*A1_coef;
    end
    for i = 1:N
        tmp = x_half(i);
        g = @(x) x - (T*(sin(x))+x-tmp)/(T*cos(x)+1);
        if X(i) <pi
            x0 = 0;
            b = New_Iter(x0,g,1e-6,1000);
        else
            x0 = 2*pi;
            b = New_Iter(x0,g,1e-15,3000);
        end
        u_true(i) = sin(b);
    end
    errors(d) = max(abs(u_true-u0));
end
orders = log2(errors(1:end-1)./errors(2:end));
%plot(x_half,u0,x_half,u_true);
plot(x_half,u0,x_half,u_true);
function [Pn, DPn] = simple_basis(x0,h,k)
Pn = @(x) ((x-x0).^k)/(h.^k);
DPn = @(x) k*((x-x0).^(k-1))/(h.^k);
end
function [u0,A0_coef] = Euler_forward(u0,A0_coef,x,x_half,k,N,h,dt,b_fun,db_fun)
New_A0_coef = zeros(k+1,N);
New_u0 = zeros(N,1);
f = @(x) 0.5*(x.^2);
for i = 1 : N
    M = zeros(k+1,k+1);
    b = zeros(k+1,1);
    alpha = 1.5;
    if i == N
        q = 1;
    else
        q = i+1;
    end
    if i == 1
        p = N;
    else
        p = i-1;
    end
    %此区间的左右端点值
    uh_left = 0;
    uh_right = 0;
    for l = 1:k+1
        Pl = b_fun{i,l};
        uh_left = uh_left+A0_coef(l,i)*Pl(x(i));
        uh_right = uh_right+A0_coef(l,i)*Pl(x(i+1));
    end
    %下一区间的左端点
    uh_right_p = 0;
    for l = 1 : k+1
        Pl = b_fun{q,l};
        uh_right_p = uh_right_p+A0_coef(l,q)*Pl(x(q));
    end
    %上一区间的右端点
    uh_left_m = 0;
    for l = 1 : k+1
        Pl = b_fun{p,l};
        uh_left_m = uh_left_m + A0_coef(l,p)*Pl(x(p+1));
    end
    flux_plus = 0.5*(f(uh_right_p)+f(uh_right)-alpha*(uh_right_p-uh_right));
    flux_minus = 0.5*(f(uh_left)+f(uh_left_m)-alpha*(uh_left-uh_left_m));
    for m = 1 : k+1
        Pm = b_fun{i,m};
        DPm = db_fun{i,m};
        b(m) = -flux_plus*Pm(x(i+1))+flux_minus*Pm(x(i));
        burger_fun = @(x) 0.*x;
        for l = 1: k+1
            Pl = b_fun{i,l};
            PmPl = @(x) Pm(x).*Pl(x);
            burger_fun = @(x) burger_fun(x)+A0_coef(l,i)*Pl(x);
            M(m,l) = quadgk(PmPl,x(i),x(i+1));
        end
        tmp = @(x) burger_fun(x).^2;
        tmpDPm = @(x) tmp(x).*DPm(x);
        b(m) = b(m) + 0.5*quadgk(tmpDPm,x(i),x(i+1));
    end
    New_A0_coef(:,i) = A0_coef(:,i) + dt*(M\b);
    tmp = 0;
    for l = 1 :k+1
        Pl = b_fun{i,l};
        tmp = tmp + New_A0_coef(l,i)*Pl(x_half(i));
    end
    New_u0(i) = tmp;
end
u0 = New_u0;
A0_coef = New_A0_coef;
end
%生成区间[a,b]上的标准Legendre多项式
function [Pn,DPn] = Legendre(n,a,b)
Q = @(x) 0.*x;
DQ = @(x) 0.*x;
R = @(x) 0.*x + 1;
DR = @(x) 0.*x;
P = @(x) ((2*0 + 1).*x.*R(x) - 0.*Q(x))/(0 + 1);
DP = @(x) ((2*0 + 1).*(x.*DR(x) + R(x)) - 0.*DQ(x))/(0 + 1);
if n == -1
    Pn = Q;
    DPn = DQ;
elseif n == 0
    Pn = R;
    DPn = DR;
elseif n == 1
    Pn = P;
    DPn = DP;
else
    for i = 3:n + 1
        k = i - 2;
        Q = R;
        DQ = DR;
        R = P;
        DR = DP;
        P = @(x) ((2*k + 1).*x.*R(x) - k.*Q(x))/(k + 1);
        DP = @(x) ((2*k + 1).*(x.*DR(x) + R(x)) - k.*DQ(x))/(k + 1);
    end
    Pn = P;
    DPn = DP;
end
%尺度变换
c = (a + b)/2;
h = (b - a)/2;
Pn = @(x) sqrt((2*n + 1)/(2*h))*Pn((x - c)/h);
DPn = @(x) sqrt((2*n + 1)/(2*h^3))*DPn((x - c)/h);
end
%生成一个cell存放所有需要的基函数
function [b_fun, db_fun] = pro_basis(k,x)
N = length(x)-1;
b_fun = cell(N,k+1);
db_fun = cell(N,k+1);
for i = 1 : N
    for l = 1 : k+1
        [Pl,DPl] = Legendre(l-1,x(i),x(i+1));
        b_fun{i,l} = Pl;
        db_fun{i,l} = DPl;
    end
end
end
function x = New_Iter(x0,g,tol,N)
x = g(x0);
n=1;
while abs(x-x0)>=tol && n<N
    x0 = x;
    x = g(x0);
    n = n+1;
end
x = g(x);
end
%k个点的高斯型求积公式
function v = Gauss_Integral(k)

end